function [matched, beta, t] = PSM(treatment, control, control_already_matched, covariates, cov_index)
    matched=[treatment,nan(size(treatment,1),2)];
    
    mkt=unique([treatment;control]);
    n=size(mkt,1);
    
    treat(1:n,1)=0;
    index=ismember(mkt, treatment);
    treat(index)=1;
    
    index=ismember(covariates(:,1),mkt);
    covariates=covariates(index,:);
    
    Y=categorical(treat);
    X=covariates(:,cov_index);
    
    % logit regression
    [beta, ~, stats] = mnrfit(X,Y);
    t=stats.t;
    
    % calculate p-score
    X=[ones(n,1),X];
    fit=exp(X*(-beta));
    pscore=fit./(fit+1);
    
    % match
    treat_index=find(ismember(mkt, treatment));
    pscore_treat=pscore(treat_index);
    [pscore_treat_new, id]=sort(pscore_treat,'descend');
    treat_index_new=treat_index(id);

    control_index=find(~ismember(mkt, treatment) & ~ismember(mkt, control_already_matched));
    pscore_control=pscore(control_index);
    
    for i = 1:size(treat_index,1)
        pscore_treat_i=pscore_treat_new(i);
        pscore_valid_control=pscore_control;

        [ctrl, ctrl_index]=min(abs(pscore_valid_control-pscore_treat_i));

        matched(id(i),2)=mkt(control_index(ctrl_index));
        matched(id(i),3)=abs(pscore_treat_i-pscore_valid_control(ctrl_index));
        control_index(ctrl_index)=[];
        pscore_control(ctrl_index)=[];
    end
    
%     temp=[mkt, ones(length(treat),1),treat];
%     temp=array2table(temp,'VariableNames',{'mkt','y','treat'});
%     temp=[temp,array2table(covariates(:,cov_index))];
%     filename = strcat('logistic', num2str(year),'.xlsx');
%     [N, T, Raw]=xlsread(filename, 1); % only clear the first sheet
%     [Raw{:, :}]=deal(NaN);
%     xlswrite(filename, Raw, 1);
%     writetable(temp,filename,'Sheet','Sheet1','WriteVariableNames',true);

    
end

